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The dispersion of inertial particles continuously emitted from a point source is analyt- 
ically investigated in the limit of small inertia. Our focus is on the evolution equation 
of the particle joint probability density function p{x,v,t), x and v being the particle 
position and velocity, respectively. For finite inertia, position and velocity variables are 
coupled, with the result that p{x,v,t) can be determined by solving a partial differen- 
tial equation in a 2(i-dimensional space, d being the physical-space dimensionality. For 
small inertia, (a;, 7;)-variables decouple and the determination oi p{x,v ,t) is reduced to 
solve a system of two standard forced advection-diffusion equations in the space variables 
X. The latter equations are derived here from first principles, i.e. from the well-known 
Lagrangian evolution equations for position and particle velocity. 
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1. Introduction 

The study of particulate matter (PM) in flowing fluids is a problem attracting much at- 
tention because of its myriad of applications in different realms of science and technology. 
In way of example, we mention the implications on: climate dynamics and hydrological cy- 
cles ( IPCC (2007)1 Lohmann & Feichter (2005)), mainly in connection to global climate 
changes originated by PM- induced cloud formation (Celani, Mazzino & Tizzi (2008)); 
environmental sciences (Shi et al. (2001)), in relation to pollution and deterioration of 
visibility; epidemiology (Dockery et al. (1993) ), in connection to adverse health effects in 



humans; and, finally, classical fluid dynamics, e.g. to understand how a flow field influ- 
ences the particle concentration ( Eaton & Fessler (1994) Balkovsky, Falkovich fc Fouxon (2001)[ 
Wilkinson fc Mehlig (2003)] |Bec (2005)[|Bec, Cencini fc Hillerbrand (2007)D . 

The complex character of this discipline simply arises from the many interactions be- 
tween totally different fields of science and technology. If, on the one hand, key achieve- 
ments have been obtained on both the epidemiological aspect and on the side of exper- 
imental studies, aiming at eliciting information on the distribution and number concen- 
tration in the environment, on the other hand very little is known about the dynamical 
aspect of particulate matter in situations of interest. As an instance, the basic equations 
ruling the spatio-temporal evolution of the particle concentration field are up to now 
unknown in the presence of paradigmatic environmental sources, such as, e.g., point- and 
line-source emissions (mimicking the release from a chimney and from a street, respec- 
tively). This seems to be the case in spite of the fact that the basic research dealing 
with inertial particles in fiuids is a field in strong development, and where important 
results have been obtained in the last few years (see, e.g. 
Ruiz, Macias fc Peters (2004)[ [Pavliotis fc Stuart (2005)[ ). 



Wilkinson fc Mehlig (2003) 



Our specific aim here is to make a first step in the intermediate framework lying 
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between the realm of abstract models of transport of inertial particles (sometimes too 
abstract to have direct implications on applicative fields) and that of empirical mod- 
els of transport (sometimes too empirical and with scarce contact with the underlying 
physics). Along this way, we will focus on the dynamics of PM when emitted from a 
point source. Note that only recently the point-source emission for the case of neutral 
particles (i.e. having the same density as the surrounding fluid) has been addressed quan- 
titatively from first principles (Celani, Martins Afonso & Mazzino (2007) I. The main dif- 
ference between the cases of neutral and of inertial particles is that, in the former case, 
the equations for the space/time evolution of particle concentration are very well known, 
and the resulting phenomenology has been elucidated only recently; in the latter case, we 
are still waiting for the governing field equations. Here, our goal is to deduce such equa- 
tions from first principles, i.e. from the known Lagrangian evolution of inertial particles 
( |Maxey fc Rile y (1983)|). 

The paper is organized as follows. In § [2] we sketch the problem under consideration 
by recalling the significant equations and by performing quantitative balances to justify 
our approach. In §[3] we analyse our small-inertia expansion order by order, focusing on 
the related equations which stem out as solvability conditions. Conclusions follow in § |31 
The appendix is devoted to briefly recalling the main steps leading to the Fokker-Planck 
equation for the phase-space density evolution, starting from the well-known Lagrangian 
equations ruling the particle evolution. 



2. Position of the problem 

Let us consider the problem of dispersion of single, small, spherical inertial particles 
emitted from a localized source, such as a chimney releasing some pollutant in the atmo- 
sphere or as an injection point (a syringe) in a microchannel, at a rate T~^. For the sake of 
simplicity, the spatial structure of the source is approximated as punctual, and located in 
the origin of our frame of reference. In terms of the release velocity, the emission distribu- 
tion (denoted by /) will be left as unspecified (with the only constraint of normalization) 
for most of the general calculations, and then for exemplificative applications will be 
modelled as a Gaussian, centered on an average value and with standard deviation a. 
The axes are chosen such as to have aligned with the positive Xd direction (our polar 
axis): for the specific case of atmospheric pollution released from a smokestack, it verti- 
cally points upwards, i.e. opposed to the gravity acceleration g, due to buoyancy effects. 
We consider as free parameters both the space dimension (even if the three-dimensional 
case should be born in mind for applicative purposes) and the coefficient /3, built from 
the particle (pp) and the fiuid (pf) densities as follows: /3 = 3/9f/(pf -I- 2pp). 
Neglecting some corrective terms but taking into account the added-mass effect in a 



simplified way (Martins Afonso (2008) I, the dynamical equations ruling the evolution of 
the particle position {X{t)) and covelocity {V{t)) in an incompressible flow u{x,t) read 
dMaxey fc Riley (1983)|): 



X{t) = V{t)+Pu{X{t),t) 

^^^^ ^ _Vit)-il-P)u[XitU ^ _ ^ V2^^^^^ (2.1) 
T r 

The Stokes time r = E? /{iv(3) expresses the typical response delay of the particle to flow 
variations, i.e. the relaxation time in Stokes' viscous drag {R and u being the particle 
radius and the fluid viscosity, respectively); ry(t) is the standard white noise associated 
to the particle Brownian diffusivity k. Notice that, when /? ^ 0, there is a discrepancy 
between particle's velocity, X{t), and covelocity, V{t) = X{t) — Pu{X{t),t); for steady 
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flows, this means that, unless pp ^ pi, the values of and tr^ should be interpreted as 

-'^lomission " fiu\x=o and X^lcmission +/?^u|^^oj respectively (here the bar represents the 
average over the velocity distribution of the particle at the emission point). 

If we consider the external (unperturbed) flow as having characteristic length scale 
and velocity L and U, the forced Fokker-Planck equation for the particle phase-space 
density p{x, v, t) reads: 



(1 - I3)u, 



(1 - P)9, 



P - 



S{x)f{v) 



T 



(2.2) 



where dt = d/dt, 9^ = d/dXf^ and = d/dv^j^. The derivation of this equation from 
(|2.ip is a lengthy but straightforward task, and is briefly recalled in the appendix. 
After performing the nondimensionalization (Martins Afonso (2008)) 



— 



U I 



T7 t 1 



t 



L/U 



we can introduce the adimensional numbers 
tU 



Stokes: St = 



Froude: Fr = 



U 



Peclet: Pe = 



LU 



and the vertical unit vector pointing downwards: G = g/ g. 

Note that the particle covelocity v has a nondimensionalization different from the fluid 
velocity u. Let us then try to explain the meaning of this nondimensionalization with 
y/lnjT, in view of our investigation in the limit of small inertia. In the phase space, this 
can be seen as resulting from a competition between the small-r limit, which would reduce 
the phase-space dimension by making v collapse onto the variety u{x,t\ and the effect 
of diffusivity, which counteracts this process. In order to have the nondimensionalized 
covelocity ~ 0(1) (i.e., a meaningful nondimensionalization) and, at the same time, the 
dimensional one ~ 0(J7) (because of the small-inertia limit, which is usually associated 
with small deviations of the particle from the underlying fluid trajectory), the following 
is a sufficient condition: 



0{v) - U 



St Pc 



0(1) 



(loosely) StPe- 0(1) 



(2.3) 



3. Expansion at small inertia 

In terms of the nondimensional variables introduced in the previous section, 
becomes: 



, ^ /Pc 



(3.1) 



-St" [At + ^u^dy\ + Sti/2 



/Pe 1 - /3 



L 



2 Fr^ 

Let us now perform a small-Stokes expansion of the particle probability density func- 



tion (PDF) and a Hermitianization of the problem in the spirit of Martins Afonso (2008) 



p(a;,-i;,i) ^ ^St"/V(a;,-y,i) , Pn{x,v,t) ^ e-'''/^Mx,v,t) . (3.2) 



Notice that, because p is normahzed to unity, eachp„ must be normahzed to dno- Plugging 
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(j3.2p into ()3.ip . we obtain a chain of equations for ipn at the different orders in St, which 
we analyze in detail in what follows. 

• O(St^i) : 

(V^ - t,2 + d) ^0 = =^ V^o = c-''"'^Ux, t) , (3.3) 

where is still unknown at this stage (we only know it must be normalized to tt^"*/^ = 
(/ Ave""" )^^), and will be found by imposing the solvability condition at 0(St"). 



+ d)'il! 



1 = e 



Co 



(3.4) 



where is still unknown at this stage (we only know it must be normalized to 0), and 
will be found by imposing the solvability condition at 0(St^^^). 

• 0{St°) : 



+2e-'''/2| [d, + (2/3 - l)u^d^ + (1 - pfVeu^] 

+v^v, [2(1 - /3)(a^u,) + 4(1 - I3)u^d, 
~2Ve-^d^d, - 2(1 - pfVeu,,u,] 
[V2Pe-i/29^ - V2(l - /3)Pei/2u^] 



(3.5) 



At this point, we have to impose the solvability condition (which would have been trivial 
at the lower orders 



j dv e~"'/2 (-y2 _ ^^2 + ^^2 = j^v V'2 (V^ - v'^ + d) e'^'/^ = 
^ = -2Aj(a;)|dt;/(7;) 
+2 [dt + (2/3 - l)u^9^ + (1 - pfVeu"] yd-ue""' 
+2 [2(1 - /3)(9^u,) + 4(1 - I3)u^d, - 2Ve-^d^,d^ - 2(1 - (ifVeu^u,] / dv v^v^c-^' 



V2Pe-^/^d,, - V2(l - 0)Pe^^\^ 6 Uvv^^e-" 



+2 

=^ {dt+u-d- Ve-^d^)£,o = A-^-d/^Six) . 
This means that, in the limit St — > 0, the nondimensional PDF is 



(3.6) 
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where behaves as a passive scalar with spatial point source; moreover, after integrating 
on the covelocity variable to obtain the physical-space density, 



P{x,t) = J dvp{x,v,t) , 



(3.7) 



in dimensional quantities one obtains the well-known equation ( Celani, Martins Afonso & Mazzino (2007) ) 



{dt+u-d- nd^)P = ^d{x) (for St ^ 0) 



Now, substituting dt£,o from (j3.6p back to (|3.5p . one gets: 



(V2-„2+d)^2 



L 



+2e-"'/2| [-2(1 - (3)u^,^^, + Ve^^d'^ + (1 - (ifVeu"] 

[2(1 - /3)(9^M.) + 4(1 - (3)u^,d. 
~2Ve-^d^d, - 2(1 - pfFeu^u,] Co 

[\/2Pe-i/'9^ - \/2(l - /3)Pei/2u^] ^i} (3 



=^ = e-'^-/2 ^) ^ ^) + (a;, t)] - 2—5{x)^{v) (3.9) 

where ^2, P2^ <yi can be derived by pushing the investigation at higher orders in St (not 
reported here), and </) satisfies the forced equation 



(V^ - w2 + d) = /(t;)c'' /2 - TT 



uV2 _ ^-rf/2„-«V2 



(3.10) 



• 0(Sti/2) 



(V2 - + d) ^3 ^ "^^^ [(^ ^ /?)Pei/2M^,5(a;)(V^ - «^)0(t;) + 2Pe-i/2^^^(^)5^5(^) 

+2e-'''/2{ {dt + (2/3 - l)u^9^ + (1 - fifVev?-\ ^ 

[2(1 - /3)(5^?/,) + 4(1 - (5)u^,^, 
-2Pe-ia^5, - 2(1 - /3)2peu^u,] 
+t;^[(...)eo + (...)6] 



(3.11) 
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The solvability condition has to be imposed again at this point: 

^ = -2^2^ j Ave-''" [(1 - P)Ve^'^u^5{x){V,, - v^)c^{v) + 2^e-^'^v^c^{v)d^5{x) 
+2 [At + (2/3 - l)u,,d,, + (1 - /3)2peu2] ^ j Ave'"" 
+2 [2(1 - I3){d^u,) + 4(1 - I3)u^d, - 2Ve-^d^,d, - 2(1 - l3 fVeu^,u,] 6 ^ Av v^v^q-"" 

+2[(...)eo + (...)e2] jAvv^e--" 
+2(. . .)^o ydt; -UpUi^wAe""' 
^ (At+u-d- Pe-i52)a = 2A/2^^-'^/2pe-i/2[a^^(3.)] y dv e~-" /\^cl){v) . (3.12) 

Let us consider two different cases, corresponding I) to the presence and II) to the absence 
of rotational symmetry in /. 

I) In this case / = /(«) that means a centered and isotropic covelocity, corresponding to 
a source emitting particles with velocity /3tt(0,t) plus isotropic fluctuations (meaningful 
for /? = if t»* = 0). In this case, then also (j) = (t){v) and the integral on the right-hand 
side (RHS) of (|3.12p vanishes because of parity, therefore = const. = and (from 



A = e-'' /'«p[\/2(l - /3)Pei/2u^ _ V2Pe''/^d^]^o ■ 

When projecting onto the physical space through (|3.7p . this term integrates to zero, and 
no explicit correction appears with respect to the tracer approximation: 

P = e-^'Co + 0(St) . 

II) On the contrary, if f{v) has a preferential direction (a mean value v^, along the positive 
Xd coordinate) one has to decompose (j) = '^i „i4'(i,m)Yi^m onto the spherical harmon- 
ics and obtains a quantum-harmonic-oscillator-like equation. As / shows degeneracy on 
the azimuthal angle, i.e. no preferential direction on the plane orthogonal to the mean 
emission, the same happens for 0, which thus reduces to a sum of Legendre polynomials 
(spherical harmonics with m — 0). With the substitution into the RHS of (|3.12p in mind, 
one can discard the isotropic sector and, if e.g. 



.1/2, 



f{v) = {27ra^)-''/\-<-^-^'^"/^''" , 



one gets: 



+ {d- l)v-'d, -l{l + d- 2)w-2 _ ^,2 + c^]^(,)(^;) = e-'/2/^^^(„) 
- iV47r(2Z + l)(2^a2)-d/2g-(i-.^)«V2-^--^/2-^^-^(_2i«w,/a2) 
(with / > 1), ji being the spherical Bessel function of the first kind. Thus, 



, I d ,9 
^'2'2+^'^ 



Kl Akk'+'-'C%+'-'\e)f^,^ik) 



Hi) 



K Akk 



l+d-l 



u 



I d 

2' 2 



+ l,k^^hi){k) 



(3.13) 
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where the values of C(i) and C(i) can be found by hiiposing regularity and normalization, 
and the constant K is given by: 

K=- 



2U 



I d 

2 



(whose value for d = 3 and / = 1 is n/2). Here, £ is the generalized Laguerre function and 

U is Kummer's confluent hypergeometric function (of the second kind) ( [Gradshteyn fc Ryzhik (1965) ). 

One can further focus on the first anisotropic sector I — 1 (the only non-vanishing when 

substituting in (|3.12p ). which implies C(i) = = C(i) for regularity^ As a result, the 

RHS of (131^ reduces to 

^]lli^^-^'''^^'Pe-'^'{d.Jix)] I dz;z;V"'/2^(i)(«) , (3.14) 
the latter integral being just a number. 



4. Conclusions 

Gathering all of our information, for the particle PDF we have: 

p{x,v,t) = e-"''^^ Uo{x,v,t) + St^^^ipi{x,v,t) \ +0{St) (4.1) 
= e-'^'^oix, t) + Sfi/^e-'^' {^1 {x, t) + «^ x 

X V2{1 - l3)Ve^/^Uf,{x, t)io{x, t) - V2Pe-^'^dM^^ *)] } + 0(St) , 

with ^0 and obeying the forced advection-diffusion equations l|3.6p and (|3.12p (simpli- 
fied through (|3.14p ) respectively. 

From (|4.ip . one can obtain information on different physical properties, or otherwise 
deduce the physical-space PDF (|3.7p directly: 

P(a;,i) -^'^/%(a;,<) + Sti/V/2^i(a;,t) + 0(St) . (4.2) 

Notice that, up to the order we have investigated (0(St^^^)), gravity plays no role. In 
other words, particle sedimentation, which is known to take place with a (bare) terminal 
velocity oc St plus flow-induced corrections, thus gives rise to a subleading effect. 

In order to give a practical example of application of our equations, let us concen- 
trate on micro-powder dispersion in microchannels. In a typical microchannel size of 
L ^ 10"'' m, inside which distilled water at temperature T ^ 320 K (density pf ~ 
988kgm~^ and kinematic viscosity v ^ h 10~^m^ s~') is flowing at the typical velocity 
U ~ 10~^ms~^, we inject lead micro-powder containing particles of size 2 10~^m 
having a density pp ^ 11.3 lO^kgm"'^. The Stokes time is r '--^ 2 10~^ s and the par- 
ticle diffusivity is k ~ ksT / {Qt^v pfR) ^ 2.4 10~'^m^s~' (from Einstein's relation, 
/cb = 1.38 10"^'^ JK~' being Boltzmann's constant). In this way, the Peclet number is 
Pe = LU/k - 4 10^ the Stokes number is St = tU / L - 2 10"^, and thus StPe = 0(1), 
according to (|2.3p . Note that, despite the fact that St is quite small, the first correction 
in our equation does appear at 0{St^^^), thus giving a more appreciable 5 10"'^ which 
perfectly falls in the present perturbative scheme. 

f Actually, (j) itself needs not be regular in the v coordinate, but the condition on the inte- 
gration constants C(i) and C(i) follows from plugging (|3.13|) into p.9|) and (|3.2|l . and from the 
consideration that p2 must be regular for v ^ Q and vanishing for v — > oo. 
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Appendix. The Fokker Planck equation associated to the Lagrangian 
description 

The derivation of the Fokker-Planck equation (|2.2p from the Lagrangian equations 
(|2.ip is briefly recalled here by exploiting the corresponding unforced (sourceless) case. 
The distribution p{x,v,t) is given by the average (on the realizations of the random 
noise r/) of the product of two Dirac deltas, expressing the probability density that, at 
time t, the particle location X and covelocity V equal the spatial coordinate x and the 
covelocity variable v, respectively: 



p{x,v,t) = (S{x ~ X{t))S{v - Vm 



(Al) 



(notice that, in this expression and in the following ones, the average acts on the random 
functions X and V, which are rj-dependent, and not on x and v). A differentiation in 
time and the use of the chain rule give: 

dp /■ dSix-Xjt)) ^m^vm 95iv-V{t)) 

m ^ \^^'^ ■ dX{t) '^^^ ~ + X{t))V{t) ■ g^(^) 

Exploiting (|2.ip and of the translational invariance of 5, but keeping the term explicitly 
depending on momentarily apart, we get: 



dp 
di 



= {{V{t)+l3u{X{t),t)} 



dS{x~X{t)) 



dx 



5{v-V{t))+5{x- X{t)) X 
d5{v-V{t))' 



dv 



ivit) 



The derivatives can now be moved out of the brackets and, inside the curly braces, the 
capital letters can be replaced by the corresponding small ones because of the deltas: 

^ - -4- (K + PM^, t)}S{x - x{t))S{v vm 



dt 



d 



(1 ^ I3)u^{x,t) 



+ il-p)gASix-Xit))div-V{t)) 



Moving the curly brackets out of the averages, applying (|A and moving to the left- 
hand side (LHS) the first two terms on the RHS, we are left with: 



d d d 
+ - — {Vf, + (3u^{x, t)) + 



dt dxu 



dvu 



{1- f3)u^ix,t) 



(1 - P)9^ 



(A 2) 



If one takes into account the source contribution, which forces the addition of a term 
S{x)f{v)/T on the RHS of (lA 2\i (this can be interpreted as the imposition of a boundary 
condition), the derivation of (|2.2p is thus complete once the stochastic term [*] is made 
explicit, which can easily be done by means of a Gaussian integration by parts (here, D 
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represents the functional derivative) (Frisch (1995)): 




dS{v-V{t))DV^{t) 
dVxit) D?/^(t) 



) 



) 



Integrating (|2.ip formally in time and invoking causality (a variation of ry at time t 
cannot cause any change in X or V at the same instant, but only at later times), in the 
last expression the former functional derivative vanishes and, in the latter, only the term 
explicitly showing r] survives: 
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